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ABSTRACT 


Optimal design of laminated composite plates with cut-outs is a topic of immense 
interest to the Aerospace industry, where the goal is to reduce structural weight. 
Several criteria have to be taken into account before the final “optimal” shape of 
the cut-out can be reached, i.e. ply failure, buckling, natural frequency, etc. Sev- 
eral studies have been reported in the literature for shape optimization problems. 
However, the final “optimal” shape will be meaningful only when constraints and 
objective functions, are evaluated with acceptable accuracy. Thus, a methodology 
to estimate error in the respective quantity of interest has to be incorporated in the 
analysis. Further, tools for mesh refinement, in order to control the error in the 
quantities of interest have to be employed in order to guarantee a reliable design. 

This study is a first step towards a reliable shape optimization package, for lami- 
nated plate structures. For a fixed higher order plate model, a simple a-posteriori 
stress recovery algorithm has been developed. The recovered stress is employed to 
estimate the error in the finite element stress, and to calculate the energy norm of 
the error in the solution. The effect of higher approximation order, and mesh refine- 
ment on the quality of the obtained solution quantities (e.g. stress components and 
displacements) is studied in detail. It is observed that higher order approximations, 
along with mesh refinement, is required to obtain pointwise stress information accu- 
rately. A detailed comparison of the stresses obtained from the recovery procedure, 
as compared to three-dimensional results has been carried out for a few sample cases. 
The recovery is seen to be accurate. The shape of the cut-out is then optimized, 
with the weight as the objective function and the first-ply failure criterion as the 
constraint. It is observed that control of the discretization error (via adaptivity) 
leads to vastly different final designs, as compared to those obtained using reason- 
ably refined meshes, but without adaptivity. It is observed that without adaptivity, 
the design obtained is unsafe, as more material removal is predicted, as compared 
to that obtained using adaptivity. 
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Chapter 1 


INTRODUCTION 


1.1 GENERAL 

“ Reliability of the constraints in shape optimization effects the 
final optimal shape. ” 

Constraint evaluation plays the main role in any optimization problem. 

If the evaluated constraints itself are inaccurate then the final optimal 
solution, obviously, represents a false picture, which is dangerous in the 
design process. Hence, before proceeding for any optimization procedure 
one should have enough confidence in the evaluated constraint. Overes- 
timation is always better than underestimation. 

Composite materials are now well established and widely used in structural applica- 
tions within the aerospace, automobile and other high performance industries where 
weight reduction and directional properties are the main criteria. Cut-outs are un- 
avoidable in structures like aircrafts. They are used for door openings and windows 
in the fuselage, flow passages, wiring connections, bolting, assembly purpose etc. 
Sometimes cut-outs are used merely to reduce the structural weight like in the web 
members of aircraft wings. Aircraft structures are subjected to complex situation of 
loadings, like Aerodynamic loads, Maneuvering loads, friction loads, etc. In general, 



these structures are weight sensitive and reliability stringent. The design of light 
weight and safe laminated composite plates with cut-outs has thus become an im- 
portant topic of research. Here, an attempt is made to optimize the cut-out shape in 
moderately thick composite laminated plates for minimum weight design with first- 
ply failure criterion as a constraint. The effect of mesh refinement on the quality of 
the computed stresses, and hence the failure criterion will also be discussed. 


1.2 REVIEW OF THE STATE OF THE ART 



Figure 1.1: Laminated plate with cut-out 

Here, an attempt is made to demonstrate the effect of reliability of constraints 
on the final optimal solution. Initially, we start with a given circular hole and 
optimize its final shape for minimum weight subjected to constraints, i.e. Tsai-Wu 
failure criterion, without considering reliability of the computed data.The process is 
then repeated with a control on the reliability of the computed data, i.e. effect of 
adaptivity on final optimal shape. 

Many researchers have attempted to optimize the composite laminate for de- 
sign variables like ply thickness and ply orientation in order to obtain minimum 
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weight designs subjected to several constraints, such as maximum strength, max- 
imum stress, first ply failure load (or reliability requirements) etc. (see [1]- [8]). 
The optimum design of laminated plates for maximum buckling load has also been 
attempted (see [9]- [11]) with constraints on the natural frequency. Others have 
worked on the optimization of composite plates with a cut-out (see [12], [13]). 
Limited literature is available on shape optimization of composite structures (see 
[14], [15]). Botkin [14] has worked on shape optimization of the stamped sheet 
metal parts with buckling and stress constraints. Sivakumar et al [15] have worked 
on optimization with dynamic constraints. Generally, the focus in all the studies 
mentioned above has been to demonstrate the effect of optimization on the final 
design. Thus, a fixed finite element mesh has been used, with quadratic approxima- 
tion, to obtain the results. The effect of the discretisation error on the final optimal 
design has not been discussed. 


1.3 LAYOUT OF THESIS 

The organisation of the thesis is given as: 

1. The first chapter deals with a discussion on the state-of-the-art and the goals 
of this study. 

2. The second chapter is devoted to a discussion of the plate model and energy 
formulation. 

3. The third chapter deals with the finite element formulation of the problem. 

4. The fourth chapter explains the recovery based a-posteriori error estimator for 
laminated plates. It is then followed by adaptive technique employed for reliable 
constraints evaluation. 

5. The fifth chapter deals with optimization problem formulation, optimization 
technique and explains the Tsai-Wu failure criterion. 

6. In the sixth chapter, first we present validation results for the plate model and 
Tsai-Wu failure criterion followed by other numerical results and inferences from 
our study. 
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7. Finally, in the seventh chapter, we have presented the conclusions made from our 
study and discussion on the future scope of this work. 
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Chapter 2 


HIGHER ORDER PLATE 
THEORY 


The analysis of laminated plates is based on the choice of a plate theory. Several 
plate theories have been developed (see [7], [8], [15]), with assumed variation of the 
displacement field in the transverse direction. Below we present the details of one 
such plate theory due to Kant et al [16]. 

2.1 DEFINITION OF DISPLACEMENT 
FIELDS 

Symmetric laminates find many applications in the aircraft industry. Although 
symmetric laminates are simple to analyse and design, some specific applications 
of laminated composites require unsymmetric laminates. For example, the coupling 
between bending and extension exhibited by this type of laminates is an essential 
feature of jet turbine fan blades with pretwist. It can be noted that the theories 
for unsymmetric laminates are applicable to symmetric laminates as a special case. 
Unlike symmetric laminates, unsymmetric laminates exhibits the membrane-flexure 
coupling phenomenon, which necessitates the use of a displacement field containing 
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both, membrane as well as flexure deformation terms which contribute to the overall 
response of a laminate. Thus, to begin with, the displacement field 

{V(x,y,z)} = [ u(x,y,z ), v{x,y,z), w(x,y,z) ] T 


is derived from the expanded Taylor’s series in terms of thickness coordinate This 
is given as (see [16] for details): 

u(x, y, z) = u 0 (x, y) + z9 x (x , y) + z 2 4> x (x, y) + z 3 ijj x (x, y); 

v(x, y, z) = v 0 (x, y) + z9 y (x, y) + z 2 <j) y (x, y) + z 3 -tp y (x, y); (2.1) 

w(x,y,z) = w 0 (x,y) 

In the expansion (2.1), it is assumed that transverse normal strain e zz is zero. 

The linear strain-displacement relationships using small deformation theory can be 
written as follows: 


£xx ~b Z 9 X ix "b Z ~b Z V’xjsi 

= v <hy + 2 Qy, y + 2 2 <l>y, y + Z 3 *l>y, y , 
ixy = U 0 , y + Z9 X ,y + Z 2 (p x ,y "f ^ , J, 

+ u 0 )s + Z @y, x + Z 2 (f> y>x + z 3 

'Yyz — 9 y + 2 z <p y + 3 z 2 %l) y + Wq :y ; 

3 Yxz — 9 X -H 2 z (j) x + 32 : 'ijj x + Wo, x 

where comma (,) denotes the partial derivatives. 


The condition that the transverse shear stresses vanish on the plate’s top and bottom 
faces (see Fig. 2.1) is equivalent to the requirement that the corresponding strains 
be zero on these surfaces, i.e. 


7j „(*,!/, ±|) = 7»(s. y, ±|) = 0 


On introduction of the conditions as given above in the expressions for transverse 
shear strains, the following relations are obtained. 


^x — ^y — 0 and 

^ y — + w 0’y) 5 V'x = ~W^ X + w 0,x) 


(2.3) 
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Figure 2.1: Representation of arbitrary three dimensional domain 

The displacement field of Eq. (2.1) is modified by setting <j) x and 4> y to be zero 
according to conditions of Eq. (2.3). The resulting displacement field is written 
below: 

u(x, y, z) = u 0 {x,y) + z 9 x {x,y) + z 3 ^ X (x,y); 

v(x, y, z) = v 0 (x, y) + z 9 y {x , y) + z 3 ip y {x, y); (2.4) 

w(x,y,z) = w 0 (x,y ) 

In Eq. (2.4) u, v and w are the displacements along x, y and z directions respectively, 
no, vo and w 0 are the mid-plane displacements while 6 X , 8 y are rotations about y and x 
axes, respectively. Whereas, ip x and ip y are higher order terms in the Taylor’s series 
expansion and are also defined at mid-plane. Thus, the generalised displacement 
vector {5} of the mid-surface contains seven degrees of freedom (DOF) and is given 
by: 


{ 5 } = {u o ,vo,wo,0 x ,e y ,'ip x ,'iljy} T . 
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The corresponding strain-displacement relationship are: 


€xx — ^(bx T Z 9 X , X + Z? Ipxix'i 
£ J W = v 0vy d* Z Qy, y + Z S 1pyi y ] 

Ixy = U 0 , y + Z 9 x , y + Z 3 7p x , y 
+V 0 , X + Z 9 y , x + 2 3 1py, x ', 
lyz = 9y + 3 2: 2 -ipy + w 0n 
7x Z = + 3 z 2 lf) x + Wo,x 

2.2 LAMINATE CONSTITUTIVE 
EQUATIONS 


A unidirectional fibre reinforced lamina is treated as an orthotropic material whose 
material symmetry planes are parallel and transverse to the fiber direction. 



Figure 2.2: Co-ordinate axes in lamina 


Following fig. (2.2) the material co-ordinates axes L and T are defined parallel 
and perpendicular to the fibre direction respectively, while global co-ordinates are 
x and y. The angle between global co-ordinate axis x and fibre direction, material 
co-ordinate axis L, is a and is known as the orientation angle. As a sign convention 
anticlockwise angle is taken as positive. 
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Generalized Hooke’s law 

In the formulation of lamina constitutive equation the following two assumptions 
are made. 

1. The lamina is a continuum. 

2. It behaves as a linearly elastic material. 

Further, at the micro-level the following assumptions are made about the material: 

1. Perfect bonding between fibres and matrix exists. 

2. Fibres are parallel and uniformly distributed throughout. 

3. The matrix is free of voids or micro cracks and initially in a stress free state. 

4. Both fibres and matrix are isotropic and obey Hook’s law. 

Stress-strain relations for the I th lamina in the material coordinate axis, whose fibers 
are oriented at an angle a with reference to the x axis is given as ( [17], [18]): 

{ Oi = [Qij]i | cj (2.6) 

where {oq} is the vector of stress components, [Qij] is the stiffness matrix, and {e 7 } 
are the engineering strain components, for the I th lamina. 

The stresses and strains in the x, y and z directions are obtained by transformation 
of the relations given in equation (2.6). The transformed stress strain relations for 
the I th lamina are given as: 

Qn Qn Qu 0 0 e x 

Q21 Q22 Q26 0 0 e y 

Qi 6 Q26 Qg 6 0 0 < 7xj , > (2-7) 

0 0 0 Q 44 Q 45 'Yyz 

0 0 0 Q 45 Q55 'Yxz . 

J 1 v ■> 1 

Using the above lamina constitutive equations and integrating the stresses over the 
laminate thickness, the stress resultants in terms of inplane forces, moments and 
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shear forces per unit length, for a laminate with NL laminae are, 


{N> 1 | {<■>} ’ 

{M} r [Ad [A] [Oil M 

• {M*} = [AF [A! [0] {«*} (2.8) 

{Q} [ [0] [0] [A] J {4 

m m 

v ✓ \ / 

(Detailed expansion of Eq. (2.8) is given in Appendix A.) 

where {eo} are the mid-plane strains, {/c} are the mid-surface curvatures and {«;*}, 
{i/}, {i/*} are higher order terms and 

^ NL m - 1 

{N } =(N x ,N y ,N xy ) = J2 (<7x,cry,a X y) 1 dz (2.9) 

i=i Jz ‘ 

NL „ ( _ 1 

{M} T = (M x , M y , M xy ) =Y1 ( Ox , Oy, a xy ) 1 z dz: (2.10) 

i=i Jz ‘ 

NL - z; _ 1 

{M*} T = (M*. M* y , M* y ) = £ / (a x , a y , a xy ) l z z dz; (2.11) 

i=i Jz ‘ 

{Q } = {Qx,Qy) I ( Txz,T yz y dz ; (2.12) 

1=1 Jz ‘ 

NL - z 

{Q*} T = (Ql, Q*y) = E (r», r yz ) 1 Z 2 dz. (2.13) 

i=i 

Here, z { are the 2 coordinates corresponding to the laminae interfaces as shown in 
fig 2.3. 

The other terms and corresponding rigidity matrices are given in Appendix A. 



Figure 2.3: Geometry of multilayered laminate 
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Thus, with the assumed displacement model, the various rigidity matrices derived 


are: 


[D m ] = Membrane; [D c ] = Membrane-flexure coupling. 

[Db] = Flexure; [D s ] = Shear. 

The set of rigidity matrices [D m ],[D c ),[Db] and [D s ] are used in forming overall 
rigidity matrix [D r ] for the laminate. 
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Chapter 3 


FINITE ELEMENT 
FORMULATION 


Triangular elements are used in the finite element approximation employed in this 
study, along with hierarchic shape functions of order p (p < 4) . The mesh generation 
is done using advancing front method based automatic mesh generator. A mesh 
generated over the plate domain is shown in fig. 3.1. Below, we give the finite 
element formulation for the chosen plate model. 



Figure 3.1: Rectangular domain with circular cut-out, meshed with an advancing 
front method based mesh generator. 
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3.1 DEFINITIONS 


Let {V} be the displacement vector defined as: 

/ \ 
u 

{V} = < v > (3.1) 

w 

K 

Let the stress and strain vectors corresponding to {V"} be {cr} and {e} which can 
be defined as: 



From the generalized Hooke’s law which relates stress components to the respective 
strain components in global coordinate system 

M = [Q] {4 (3.3) 

where material stiffness matrix [Q] for orthotropic material is as given in Eq. (2.7) 
for each lamina. 

For any elastic body, the linear strain-displacement relationships using small defor- 
mation theory is given in Eq. (2.5). 

The above relations can be expressed in matrix form as: 

{e} = [D] {V} (3.4) 
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The components of displacement can be written in terms of the seven unknown 
in-plane functions, which can be written as: 

{ 5 } t = {uo,vo,wo,6 x ,e y ,i/j x ,ipy} 


The displacement vector is written in matrix form as 


{V} = [$] {{} 


(3.6) 


where 


1 0 0 z 0 z 3 0 

[$] = 0 1 0 0 z 0 z 3 

0 0 1 0 0 0 0 


(3.7) 


In the finite element approximation the functions Uo, v o, ... etc are approximated 
using n shape functions per element. 



W = £ M {^} (3.9) 

1=1 


where, n — (p + 1) (p + 2)/2 is the number of independent coefficients in an ele- 
ment, p is the approximation order, Ni is the i th shape function associated with i th 
independent coefficient in terms of normalised coordinates £ and rj and {5j} is the 
generalized displacement vector corresponding to i th independent coefficient of an 
element. 


Remark: When Lagrangian shape functions are used, the independent coefficients 
correspond to physical nodes. 
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Hence displacement vector can be written as 


{V} = [$] [AT] {d} (3.10) 

where {d} T = {{di} 7 , {^} T , and [N] is a matrix in terms of n inplane 

hierarchic shape functions given in Appendix B. 

Hence, the strain vector can be written as 

{4 = ([51 m [IV]) {d} (3.11) 

and the stress vector 

{a} = [0] ([^] [*] m {d> (3.i2) 

3.2 FINITE ELEMENT FORMULATION USING 
ENERGY PRINCIPLE 

The total potential, for the plate is given by 

n («) = xy w, (3.13) 

e=l 

where 7 r e is the total potential of the non-intersecting (but adjacent) sub-domains 
e which are part of the domain (N sub-domains are considered here). The total 
potential can be expressed in terms of internal strain energy and external work 
done W( e \ as follows: 

7T e (5) = U {e) - W {e) . (3.14) 

Strain energy of the laminate can expressed as follows: 

UM = \Iv M W r H<^ W . (3-15) 

Work done by the applied external transverse load is, 

= (3 . 16) 

= Ir+m Wo f + dA (e) -1- Sr-m w 0 f- dA (e) 
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where / * is the transverse load on the top face R + and f~ is the transverse load on 

the top fare R~ . 

1 Ik* exact solution {!’„} to this problem is the minimizer of the total potential II. 

This ran be obtained as: 


8 ^ n = o, 


whirls is also the Virtual Work Formulation of the problem in terms of the com- 
ponents of S. From this, the seven coupled equilibrium equations in terms of the 
components of 5 can be obtained. This leads to the generalised finite element for- 
mulation: 

[K] {d} = {F} (3.17) 

where (A'j is the global stiffness matrix and {F} is global load vector. 

In the next section, we are going to employ the virtual work formulation to derive 
the finite element formulation of this problem. 


3.3 COMPUTATION OF ELEMENT STIFFNESS 

MATRIX 


The stiffness matrix corresponding to assumed deformation state of an element can 
be defined by expressing the internal strain energy in terms of unknown nodal 
displacements. In the formulation of unsymmetric laminates the membrane, the 
flexure, membrane-flexure coupling and shear strains contribute to strain energy. 
The set of Eq. (2.5) can be used along with Eq. (3.12), (3.15) to express these 
strains in terms of nodal displacements. 

The internal strain energy of an element (given by area A e ) can be determined as. 

U e = | IaU) ( i d ) T ( [An] [Bra] ) M + {d} T ( [B b f [D b ] [B b ] ) {d} 

+ {d} T ( [Bsf [Ds] [Bs] ) {d} + {d} T ( [Bm] T [Dc] [B b ] ) {d} (3-18) 

+ {d} T ( [B b ] T [D c ] [B m ] ) {d} ) dA 
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Heir the first trim on the right-hand side of Eq. (3.18) is the in-plane contribution, 
the second and third terms are out-of-piane bending and through-thickness shearing 
contributions, respectively. The last two terms are contributions from coupling 

between in-plane and out-of-plane actions. 

The -strain energy expression can be written in a concise form as: 


P e) = i {<f } r {K'\ {d‘ } (3.19) 

where |/iT e ] is the stiffness matrix and {d e } are the coefficients corresponding to 
the finite element solution in an element e. All the terms in equation (3.18) are 

evaluated individually and then summed to yield [K e \. 

It is given by: 

P"] = Sa- [ ( PmF Pm] Pm] ) + ( W [A] Pd ) + ( PJ T Pd Pd ) 

+ ( P„F PJ Pd ) + ( Pd T Pd Pm] ) ] dA 

Numerical integration is carried out to get the element stiffness matrix. 

The matrices [B m ], [B fc ] and [B c ] and vector {d} are given in Appendix B. 


3.4 COMPUTATION OF ELEMENT LOAD 
VECTOR 

For the extension and/or bending problems of laminated plates, the applied external 
forces have been considered to be of the following types: 

1. Uniformly distributed load acting over the element in the ^-direction on top or 
bottom bounding planes of the plate. 

2. Sinusoidal distributed load acting over the element in the ^-direction on top or 
bottom bounding planes of the plate. 

3. TVaction load acting along the edges in the x or y-direction . 

The total external work done by these forces on an element can be expressed as 
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follows: 


= M T j s (W [Jv]f{r,} dS 


( 3 . 21 ) 


or 


W c = {(if {F e } 


( 3 . 22 ) 


where {F f } is the element load vector, and the traction vector is given as: 


{Tr} = 


T x 

T 

x y 


> 


T z 


( 3 . 23 ) 


When the loading is in-plane, 


{Tr} = < 

s. 


T x 

Ty 

0 


> 


When the loading is in the transverse direction, 


{Tr} = 




0 

0 

T z 


> 


Transverse loading, T g , can be sinusoidal as T z = P c sin( 2 —^) sin( 2 ^p) or uniform 
pressure T x = P 0 . T s , T y can be uniform inplane tractions along x and y, respectively. 


3.5 GEOMETRIC APPROXIMATION 

3.5.1 STRAIGHT EDGE ELEMENTS 


The geometry is expressed in terms of the shape functions as: 


x 

y 


1=1 


Xi 

Vi 


( 3 . 24 ) 
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wlu'iv, < V are th<* coordinates of the i th node of the element, and Ni are the 

l ^ J 

linear shape iunctions of the element, 

LINEAR MAPPING: 



Figure 3.2: Linear mapping in two dimension 

Linear mapping in general is used when all sides of mapped elements are straight 
lines (see fig. 3.2). Consider the mapping of variables from x,y to £, r] such that 
{.r. y G A} on substitution of Ni these can be written in the matrix form as 


Clearly, 


X — X\ 

y-y i 


2 2 “ 2 i 
1/2 “ 1/1 


2 3 - X\ 

2/3 Vi 



^=X 2 ~Xi ^ = 23 - XI 

ff = 2/2 ~ yi % = Vz-V\ 

On inverting the matrix it is easy to see that 

( \ 
X “ Xi 
< > 

y-yi j 


jc 

1 

{yz - yi) -(23 - 21) 

bJ 

'“A 

— (y2 — 2/1) (22-21) _ 


where, 

Hence, 


A = (X2 - ®i) (2/3 - 2/1 ) - (23 - 21) (2/2 - 2/1) 


dx A 

§r _ j/2-yi 
dx A 


Hi — g 3 -ai 

dy ~ A 

dji ai 

dy ~ A 


(3.25) 


(3.26) 


(3.27) 


(3.28) 
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Transforming d.rdy into d£,diy. 

L(>t ,h< ’ ,llfh ' r< ’ u, i ai area dxd V is formed through vectors dx and dy with magnitude 

d. \ and direction normal to the elemental area is k 

dA = dx dy = [dx x dy]-k (3.29) 


dt <ly = [ ( f? i + | <m j ) X ( | d( i + |2 dr, j ) ]. k 

= ]«*) ( 3 . 30 ) 

= |J| d£ dr) 

or the Jacobian matrix is given by 

3x dy 

j = dt 

dx dy 

m dy dy 

Substituting from equation (3.5) 

j J| = (x 2 - Xi) (t/ 3 - y{) - (a* - a?i) (y 2 - yi) 


( 3 . 31 ) 


( 3 . 32 ) 


Hence f t , f (x, y) dx dy may be written as 

J e f (*i y) dx dy = j* f (£, 77) d£ dr) 


( 3 . 33 ) 


such that / (£,v)d£dr)r=\J\f (®(f, 77), y(f, 77)). 

Numerical Integration: 

To compute f Am f (£, 77) d£ dr) Gaussian quadrature for two-dimentional integrals 

is used through which 



/(&»?) d ZdT)&Y^ f(€i,Vi)wi |J 1 

i=i 


( 3 . 34 ) 


where, n is the number of quadrature points, is the weight at quadrature 
point (&, r)i) and the Jacobian \J\ is the area of the element e. 

Order of numerical integration is greater than or equal to (p + l)/2 where, p is 
the order of the integrand. 
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Shape Function Derivatives: 


Shape functions are defined in the master element i.e. Its deriva- 

tives with respect to x and y can be given as 


Ni 


*.* 


<?a\ ag , dN, dn 

d£ dx ’ dr] dx 


Ni 


iV 


BNj d£ , dN , By 
d£ By di) dy 


(3.35) 


Since shape functions are defined in terms of normalised coordinates £ and rj the 
terms and ^ can be evaluated at any point in the domain. 


3.5.2 CURVED ELEMENTS 

The cut-out profile is represented in terms of a curved edge. Thus, at a cut-out 
edge, linear mapping will be inadequate unless a very fine mesh is employed at the 
cut-out. Further, for higher order approximations (i.e. p > 2), significant errors in 
the solution quantities can be obtained for the elements at the cut-out, if proper 
geometric transformations are not employed. Special care has to be taken at the 
curved boundaries, to ensure that the errors due to geometric approximations do 
not become significant. 


3.5.3 MAPPING CURVED ELEMENTS 

For elements with curved edges, a quadratic blending function type mapping is 
employed. The mapping can be given as: 

X = X\ Ni + X2 Ni + xz Nz + cx. Ni+ 3 ^ 

Y = yi Ni + 2/2 N2 + 2/3 N3 + N hz 

where a, J3 are suitable coefficients; * = 1, 2, 3 is the curved edge. Below, we describe 
how a and /? are obtained. 

Let us, for example, consider the curved edge to be a segment of an ellipse. For any 
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Y * 


1 



X 

Physical element 



Figure 3.3: Curved element mapping 


point on the ellipse we have the relations 


x* — a cos 9 
y* = b sin 9 


(3.37) 


The curve between two nodes is approximated by using the quadratic shape func- 
tions. Tims, we have 


x = x 1 N l (t) + x 2 N 2 (Z) +ajV 3 (£) 
V = yi NxiO + y 2 N 2 (0 +PN 3 (0 

where 


(3.38) 


iV t (£) = N 2 (£) = Nz(£) = — v/6 (1 — £ 2 ) (3.39) 

are the one-dimensional hierarchic shape functions defined on the curved edge. 

Let us further define the following functional on the edge: 

J = J l ^x-x*) 2 <% (3.40) 

or is the minimiser of J. Thus a can be obtained from 




Hence, 

a f JV,(f ) 2 = J\x'-x)N 3 dt; (3.41) 
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(3.42) 


- - J-, ( *• - * ) AT, g 

/I. «?<«{ 

Simihirly. we can And fl, by defining J = f‘ t ( y - y - ) ^ and /? as minimise!. 

3.6 BOUNDARY CONDITIONS 

1 loin the variational formulation of the problem several combinations of the ap- 
plicable boundary conditions in terms of forces/displacements, moments/rotations, 
higher moment, s/curvatures are possible. However, for simplicity, we list certain 
typical boundary conditions (considered in this study) below. 

1. Point support boundary condition 

2 . Simply supported boundary condition 

3 . Clamped edge boundary condition 

4. Symmetry edge boundary condition 

5. Free edge (all traction components are zero). 

While applying boundary conditions we look at those degrees of freedom which are 
to be set zero (see Table (3.1) and Table (3.2)). 


Point Supported 

Simply Supported 

Clamped Edge 

Symmetry Line 

Wq = 0 

Wq -0 

Uq — 0 X — ip x — 0 

vjo = 0 

U'O = 0 X = ipx ~ b 

Vo = 6 y = lj) y - 0 

Vo = 0 

6 y = Vv = 0 


Table 3.1: Typical boundary conditions on x edge ( y = constant) 


Point Supported 

Simply Supported 

Clamped Edge 

Symmetry Line 

Wo = 0 

w 0 = 0 

Vo = 6y = 'lpy=0 

wo = 0 

uq = 9 x = ip x -0 

vo = O y = ipy = 0 

£ 

II II 

o 

a 

ii 

o 


Table 3.2: Typical boundary conditions on y edge (x = constant) 
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Solving t ho finite ('lenient problem for the seven DOFS, we can obtain the ap- 
proximate values ot displacements, strains and stresses at any point in a lamina. 
However, a-priori we have no knowledge of the accuracy of the solution quantities. 
A good approximation of the stress components is essential to accurately find the 
first -ply failure criterion. In the next section, we outline a simple procedure for the 
a-posteriori recovery of a “better” stress field from the finite element solution. 
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Chapter 4 


A-POSTERIORI RECOVERY OF 
POINTWISE STRESSES AND 
ERROR ESTIMATION 

4.1 INTRODUCTION 

In a typical engineering analysis a mathematical model for physical problem is first 
selected, such that it incorporates the essential features of the actual physical prob- 
lem. The finite element method determines an approximation to the exact solution 
of the mathematecal model. The computed solution should be compared with exact 
solution of the mathematical model which is being solved. Hence, the computed 
results can be used to make engineering decisions only when one can guarantee that 
finite element solution is sufficiently close to the exact solution of the mathematical 
model problem. In general, the results of interest obtained from the finite element 
solution can be very different from those corresponding to the exact solution of the 
mathematical formulation and can lead to serious design errors. Thus, it is imper- 
ative to accompany any computational analysis with an accurate indication of the 
error in the quantity of interest. 
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The error estimator should be reliable. Unreliable error estimates are dangerous 
bemuse they could lead to a misleading confidence in the computed quantities. The 
reliability of the error estimator has to be understood with respect to the solution 
quantity of interest. Here, the error estimator is constructed for the energy norm of 
the error and it is assumed that if the error in the global energy norm is low then 
all solution quantities of interest are also reasonably accurate. 

In tills study a simple procedure for the recovery of strains, from the finite element 
solution, using patchwise data is proposed ( as an extension of the method in [20]). 
This recovered strain will be used to design a simple error estimator. 


4.2 PROCEDURE 


The idea behind the error estimation of laminated plate is that, the in-plane strains 
1 1 he recovered one) are assumed to be of order p whereas the transverse shear strains 
are an order higher i.e. (p + 1) and they are of the same mathematical form as that 
for the exact solution of the model as in Bq. (2.5) and are given as: 


(4.1) 


£ XX 

_ ,*o 

- 6 XX 

+ z e* l xx + 


f * 

t w 

= e*° 
c yy 

+ * + 2:3 

c yy 

7*xy 

= 7*°j , + * 7* xx + 2 

3 

I xy 

7 V 

— ry*^ -1_ y} 

~~ 7 yz^ * 1 yz 


7*XJZ 

— 7 xz ' 16 1 xz 



In-plane strains of the strain tensor are approximated by p th order polynomial, 


(p+J)(p+*)/* 

e *xx — e Qi 

i=l 


(4.2) 


n * *1 t ^*2 ^*0 ^*2 f *0 f ~~ 

Similarly, € xx* € xx> e yy > 6 yy> e yy' e xyjt xy 


.*2 


strain components can be approxi- 


mated as above. 


And, transverse shear strains of the strain tensor are approximated by (p+1)" 1 order 
polynomial. 
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(4.3) 



(j>+2)(ji-fS)/2 


E 

i=J 


^*G 

^ yz,i Qi 


Similarly, y*”.. strain components can be approximated as above. 

where q, art* Tiumomials (defined in terms of a local coordinate system described 

later). 



A patch is constructed by taking element J and one layer of elements around it, as 
shown in fig. (4.1). Let the centroid of element J is (X C ,Y C ). 
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A lurai nionlinalt' systim can be defined with X C ,Y, as the centre as: 



[J 

v c 1 


UJ 

1 | y- 

n J 

r 

are given as, 




<I\ ~ 1 

(12-X 


93 = y 

A = x 2 

Qs = xy 


96 = f 

97 = x 3 

aO a 

q& = ary 


99 = xy 2 

CO 

II 

O 

5 

9n = x 4 


912 = X z y 

9is = x 2 y 2 

Qu = xy 3 


9is = y 4 

9ic = x n 

9i 7 = x 4 y 


918 = X 3 y‘ 

9i9 = x 2 y 3 

920 = xy 4 


921 = y 5 


(4.4) 


(4.5) 


Now to recover a smoothened strain field we should find smoothened strain compo- 
nents, that is, 9 [(p+l)(p+2)/2] + 4 [(p+2) (p+3)/2] unknown coefficients. To get 
these coefficients, as aforesaid, the strain recovery procedure uses the principle of 
minimization of energy norm of error i.e. the energy due to errors in strain and 
stress components, over the patch considered. A typical patch over element J is 
shown in fig. (4.1). In this, the strain components of the finite element solution 
are known. The material properties and other relavent information about patch is 
available. From this, the strain energy of the error can be computed as, 

*/ = x f f~ h ( &fe ~ £* ) ' Q ( IfjE — §* ) dz dA (4.6) 

* J A pa tch 

where §* and gpp are recovered and finite element strain vectors and Q is material 

stiffness matrix, respectively. 

The minimization of this function with respect to each unknown coefficient of re- 
covered strain term gives 9 [(p+l)(p+2)/2] + 4 [(p+2)(p+3)/2] equations. Solution 
of this system of equation gives the unknown coefficients. Using these coefficients 
smoothened strain components are constructed over the patch. 

The above step is followed by Re- construction. In this, the in-plane stresses are 
constructed using material properties i.e. from constitutive equation. Although, 
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f 1» h < slie.u .strains are improved over the finite element strains, equiiiDrium 

i ( [uatiun> an usad to get the transverse stress field as the model will not give good 
transverse stress components (because strain continuity has been assumed). 

The strain energy from recovered and finite element strains are calculated over the 

whole domain, 

t'recov ~ | EjLl fj 2.* ' £* dVpatch 

rr , „ (4-7) 

t fem — 5 fj 2/ e ' £/e dV patc h, 

wlitne f and L/ em are the strain energies over the whole domain from recovered 
strains and finite element strains, respectively, a* and <r /e are recovered and finite 
element st resses, c* and £/ e are the recovered and finite element strains and N is the 

number of elements in the domain. 

The energy norm c of the error is calculated as: 


e = i feu(u* -M/ e ) 


(4.8) 


4.3 ADAPTIVITY 

I'he energy norm of the error obtained using the recovered stress field is used to 
refine the mesh. The procedure involves the computation of the energy norm of 
the error in each element, followed by a ranking of the elements in the order of 
highest contributions to the total error. The elements contributing eighty percent 
of the maximum error are refined. This procedure is repeated till convergence to 
within the specified tolerances is obtained. It should be noted that adaptive analysis 
requires repeated solution of the boundary value problem, with the modified meshes. 
Thus, the cost of computation increases (see [20] for details). 
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Chapter 5 


DESCRIPTION OF 
OPTIMIZATION PROBLEM 


The plate with dimensions of X, Y and a centrally located cut-out of elliptical shape 
with initial size 2a, 26 and its rotation is shown in fig. (5.1). Here, 2a is major 

axis, 26 is minor axis. 



x 

Figure 5.1: Initial shape of plate and cut-out 
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1 he objective of the optimization probem is to minimize the weight of the plate. 
Since all the laminae of same material are taken, the objective becomes the mini- 
mization of the plate material area. Hence, the design parameters are, semi-major 

axis - n, semi-minor axis - b. 

Objective function: Minimize W 

( )r Minimize Area of the material A = XY — nab 

Subject to Stress state < 0.8(Tsai Wu failure criterion). 

| < a,b < 04 Y 


5.1 OPTIMIZATION TECHNIQUE 

Complex Search Method algorithm is used for solving this optimization problem 

(21]. 

The algorithm begins with a number of feasible points created at random. If a point 
is found to be infeasible, a new point is created using previously-generated points. 
Usually, the infeasible point is pushed towards the centroid of the previously- found 
feasible points. Once a set of feasible points is found, the worst point is reflected 
about the centroid of rest of the points to find a new point. Depending on the 
feasibility and function value of the new point, the point is further modified or 
accepted. If the new point falls outside the variable boundaries, the point is modified 
to fall on the violated boundary. If new point is infeasible, the point is retracted 
towards the feasible points. The worst point in the simplex is replaced by this new 
feasible point and algorithm continues for next iteration. 

Here, the reflection parameter used is 1.3 and convergence parameter is 0.01. 
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5.2 FAILURE CRITERIA 


5.2.1 TSAI-WU FAILURE CRITERION 

l ht* failure criterion based on interactive failure theories is considered. We have 
used here second-order tensor polynomial criterion proposed by Tsai-Wu [22]. It is 
a complete quadratic tensor polynomial with the linear terms included. The most 
compact form for expressing this theory is through tensor notation: 


Fi c n + F (j &i a 3 < 0.8 t, j = 1, 2, ....6 


(5.1) 


where F, and F^ are the strengths tensors established through experimental proce- 
dures and are related to failure strengths in principal lamina directions. Here, right 
hand side of above equation is deliberately taken 0.8 instead of 1 as a safe constraint 
for the optimization problem mentioned earlier. For orthotropic lamina this reduces 

to: (see [27], [28]) 


f'lO'i + F2O2 + FgOg + F n al + F22C2 + F33O3 + F 44 <t| 

+ ^ 55^5 + FeeOg + 2F\20 r i<J2 + 2Fi§oi<j§ + 2F2QO'2C r 6 5: 0.8 


(5.2) 


and for plane stress this reduces to: 


Fl&l 4 ^2^2 + FqCTq 4 Fllo\ 4 F 22&2 d* FqqCq 
4*21^12^1^2 + 2jF\$<T i<Jg 4 2jp26^2^6 ^ 0*8 


(5.3) 


The strength tensors are given below: 

= F 2 =4 _ 4 F% — ~5c’ > 

= XtXc ’ ^22 = Y T Y C ’ ^ 33 — Z T Zc> ( 54 ) 

F44 ” F55 = Fgg 

F n ss -I ^4^’ i?13 = “2 ^ = “2 vu T ycZ T z c 

where subscript T denotes tensile strengths, C denotes compressive strengths. X,Y 
and Z are the strengths in L,T and T' directions, respectively. R,S,T are shear 
strengths in TT\ LT? and LT planes, respectively. T' is pependicular to plane LT. 

(see fig. 2.2). 
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5.2.2 PROCEDURE FOR FIRST-PLY FAILURE 


ANALYSIS 

The first-ply failure analysis is based on the assumption that a given ply would fail 
if t he failure index at any point within the ply reaches a value of unity. The failure 
index is defined as ‘F, cri + F,j a ,■ aj\ The steps used in the iterative procedure are 

listed here: 

Step 1 Find the displacement field for a small initial load. 

Step 2 Find the stresses for the above displacement field. 

Step 3 Transform the stresses into the material coordinates. 

Step 4 Find the maximum failure index. 

Step 5 Check whether the laminate has failed or not. 

Step 6 If not, increase or decrease the initial load appropriately. 

Step 7 Repeat the steps 1 to 6 until the laminate fails. 

The maximum failure index is determined by carrying out a sequential search at 
certain preselected points within the laminate. The sequence of search is as follows: 

1. Consider the first element of the finite element mesh. 

2. Consider the first Gaussian point within the element. 

3. Consider the first ply (from top) of the laminate. 

4. Consider the top of the current ply. 

5. Find the failure index at that location. 

6. Check whether this is greater than the previous value. 

7. If yes, store the element number, Gaussian point number, ply number, 
location on the ply and failure index and go to the next location (middle of the 

ply, bottom of the ply, etc). 

8. If answer in Step 6 is no, then do not store the failure data, but go to the next 
location and repeat Steps 5-8. 

9. Go to the next ply. 
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Go to I lit* next Gaussian point. 

Go to next element. 

Continue the search until all the elements in the finite element mesh are 

searched for t he maximum failure index. 
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Chapter 6 


NUMERICAL RESULTS AND 
DISCUSSION 


The shape optimization process requires computation of constraint information (e.g. 
first-ply failure criterion). In order to guarantee the reliability of the constraint infor- 
mation, the accuracy of the plate model (and its implementation), the discretisation 
procedure and postprocessing module have to be established. In this chapter, a 
detailed validation study for the plate model, along with the influence of mesh re- 
finement and p - inrichment, is done. This is followed by the shape optimization 
problem with and without adaptivity. 


6.1 VALIDATION OF PLATE MODEL 

In this section, first we define the problem for validation study. It is followed by 
validation study for Thin, Moderately thick, and Thick plates. Lastly, validation 
for first-ply failure load of laminates is presented. 

Remark 1: The definition thin, moderately thick and thick plates depends upon 
length to thickness (f ) ratio, loadings, boundary conditions etc. Here, we have 

defined them roughly. 
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!h mark Unless specified, the analysis is done for full plate. 


6.1. 1 PROBLEM DEFINITION 


yf 

Y 


I 


x 


jL 

12 


h 




T ; 

x 

- 



i st . piy 

2 nd ply 
N th ply 


Figure 6.1: Problem plate for plate model validation 

A laminated plate of the dimensions shown in fig. (6.1) has been considered for 

analysis. 

The plate is loaded transversely on the upper surface through the sinusoidal load: 

. . . .m7rx . ,n-K y. / fil N 

— 9c sin( ) sin( y ) v • i 

or with uniformly distributed transverse load: 

<iz{ x >y) — ~Qo ( 6 - 2 ) 

The transverse mid-plane displacement at the centre of the plate, w 0 (unless speci- 
fied) is normalised as w 0 : 

w 0 - wo ff^r 10 2 

When load is uniformly distributed transvese load, the term g e is replaced by g 0 . 
Unless specified in all problems the mesh types used are shown in fig. (6.2). 
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MESH II 



MESH III MESH IV 


Figure 6.2: Mesh patterns used for plate validation study 

6.1.2 THICK PLATES 

The plates with X/h ratio upto 8 are treated as thick plates. 

PROBLEM 1 

A three layer simply supported square laminated plate under sinusoidal transverse 

load is analysed. The lamina properties are as follows: 

= 25; E 2 = E z = 6.89 x 10 3 N/mm 2 -, 
gaa — n q- = Qp- = 0.5; 

E* ~ U Z> Ei E 2 u 

= ^23 =: ^13 == 0.25; 

X = 25.4 mm; f = 4; 

Thickness of each layer = | ; Fibre orientations = [0°/90 /0 ]; 
q c = 6.89 x 10 ~ z N / mm 2 ] m = 1 ; n = 1 

Results of this problem are tabulated in table 6.1, 6.2. 
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Mesh Type 

Wo 

ri p x io - 5 

Mesh I 



-1.9230 

-1.7371 

Mesh II 

-1.9242 

-1.7395 

I Mesh III 

-1.9258 

-1.7397 

Mesh IV 

-1.9260 

-1.7397 

Kant et al [16] 

-1.9058 

- 


Table (5.1: Validation study for thick plate with mesh refinement and p = 3 


V 

Wq 

rip x io -5 

1 

-1.8339 

-1.6009 

2 

-1.9267 

-1.7388 

3 

-1.9258 

-1.7397 

4 

-1.9261 

-1.7397 

Kant et al [16] 

-1.9058 

- 


Table 6.2: Validation study for thick plate with p refinement for Mesh III 

We ran see from above results ( see table (6.1), (6.2) ) that: 

• The results from our study are close to the one with reference results osbtained 

by Kant et al [16]. 

• For p= 3 the effect of mesh refinement shows that results converge with mesh 

refinement above III. 

• And with p refinement it converges for p > 3. 

• The reference results are for a fixed rectangular mesh with four elements and 
with no convergence analysis. Hence, we see difference in the the two results. 
The total potential converges to a value of -1.7397 x IO -5 . 
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PROBLEM 2 


A i* - ( tangului laminated plate simply supported along all edges and subjected to 
n unsxerst' vmusutdal loading is analysed. The laminae properties are listed as fol- 

i< tWM 

Ei 138 N/mm 2 ; E 2 - 9.3 N/mm 2 ; 

= ("; 21 = 4.6 N/mm' 2 : G u ~ 3.1 N/mm 2 ] 
i'vji --- 0.3; i/j3 — i/23 = 0.5; 

,Y — 5 mm: V = 3.5 mm; /i = 1 mm; 

Thickness of each layer = Lamination scheme [0/90] s ; 

«/ r 1 N/mm 2 ; m = 1; n = 1 

The results are analysed for quarter plate with mesh shown in fig. (6.3) and are 

t abulated in table 6.3 through 6.6. 


Mesh Type 

w 0 x 10~ 4 

n p x io- 4 

Mesh I 

-5.8090 

-3.1683 

Mesh II 

-5.8439 

-3.1954 

Mesh III 

-5.8467 

-3.1974 

Mesh IV 

-5.8468 

-3.1975 

3-D solid [19] 

-5.7363 

-3.0851 


Table 6.3: Validation study for thick plate with mesh refinement for p = 2 


V 

wo x 10 4 

lip x 10" 4 

1 

-5.2480 

-2.8603 

2 

-5.8439 

-3.1954 

3 

-5.8467 

-3.1975 

4 

-5.8469 

-3.1975 

3-D solid [19] 

-5.7363 

-3.0851 


Table 6.4: Validation study for thick plate with p refinement for Mesh II 
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Mfdi 

It 

<t*y 

1 xy 


a v 

T xy 

Type 

i f;, - h) 

f it b h\ 

V ' 2* 2 ' 

M,-f) 

/a b h\ 

V 2 ’ 2 1 2 * 

fab h\ 

V 2 ’ 2 ’ 2 / 

(a, b, -|) 

i. i ; 

1 

1.8848 

-0.9662 

8.0526 

1.7647 

-0.8916 

Mi . ii 11 

7.9917 i 

1.7597 

-0.9057 

7.9523 

1.7496 

-0.8990 

M.-h HI 

i I 

7.9137 

1.7410 

-0.8989 

7.9118 

1.7403 

-0.8982 

Mesh IV j 

7.9042 

1.7384 

-0.8980 

7.9042 

1.7382 

-0.8970 

:M> solid ;i!V i 

■ ■ i 

7.7388 

1.9267 

-0.8602 

7.7388 

1.9267 

-0.8602 

Table 6 

.5: Comparison of stresses with mesh refinement for p 

= 2 

i v | 


a y 

r\y 

Vx 

j 

r xy 


/n h h \ 

*2 T 2' 2 ' 

b _ A\ 

v 2 ’ 2 T 2 * 

f) 

( a b h\ 

v 2 } 2 J 2 / 

(a b h\ 

V 2 > 2 5 2 ' 


! 1 

7.1042 

1.5524 

-0.8073 

6.9450 

1.5100 

-0.7813 

2 

7.9947 

1.7597 

-0.9057 

7.9523 

1.7496 

-0.8990 

| 3 

7.8687 

1.7298 

-0.8936 

7.8972 

1.7363 

-0.8970 

4 

s 

7.8663 

1.7513 

-0.8980 

7.8980 

1.7365 

-0.8972 

| 3-1) solid [19] 

7.7388 

1.9267 

-0.8602 

7.7388 

1.9267 

-0.8602 


Table 6.6: Comparison of stresses with p refinement for Mesh II 


In the above tables 6.5, 6.6 a*,a*,r* y denotes recovered stresses and a x ,a y ,r xy 

denotes finite element stresses. 

From above results we observe that: 

* The results converge for mesh refinement above III. 

* With p refinement for mesh III results converge for p> 3. 

* The plate model considered shows the total potential values higher than the 
reference values. This difference between total potential values is because 
the plate model considered here assumes a smooth strain field through the 
thickness, which leads to a jump in the values of the transverse stresses at 
the interfaces. This is in direct violation of the stress continuity requirement 
(from 3-D solid elasticity) at the interface. 
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MESH III 



MESH II 


a 

69 

a 

89 

a 

sia 

Si 

a 

69 

ai 


m 

fiS 

IB 

as 

IB 

asiB 

at 

IB 

at 

IB 

at 

m 

a* 

IB 

as 

IB 

asiB 

at 

IB 

at 

IB! 

101 

IB 

as 

IB 

an 

IB 

atiB 

at 

IB 

as 

IB 

at 

a 

* 

a 

Cl 

a 

*a 

89 

a 


a 

89 

tea 

as 

PB 

as 

IB 

atiB 

ai 

IB 

as 

IB 

at 

a 

69 

a 

fin 

a 

ssa 

89 

a 

K 

a 

B 

ib 

as 

IB 

as 

IB 

asiB 

at 

IB 

as 

IB 

at 

IB 

an 

IB 

as 

IB 

asiB 

at 

IB 

as 

IB 

as 

a 

* 

a 

* 

a 

ssa 

IK 

a 

89 

a 

69 


at 

IB 

at 

IB 

atiB 

at 

IB 

at 

IB 

at 

[ib 

at 

IB 

as 

IB 

at ib 

as 

IB 

as 

IIB 

as 


MESH IV 


Figure 6*3: Mesh patterns used for plate validation study 
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6.1.3 MODERATELY THICK PLATES 


Hii’ plates with X f h ratio between 8 to 15 are treated as moderately thick plates. 

PROBLEM 3 

The problem description is same as problem 1 with the change y = 10. 

Hesults are tabulated in table 6.7, 6.8. 


Mesh Type 

^0 

n p x 10“ 4 

Mesh I 

-0.7174 

-1.0098 

Mesh II 

-0.7171 

-1.0127 

Mesh III 

-0.7176 

-1.0129 

Mesh IV 

-0.7176 

-1.0129 

Kant et al [16] 

-0.7164 

' 


Table 6.7: Validation study for moderately thick plate with mesh refinement for p 

3 


P 

Wo 

a 

X 

i — L 

o 

i 

1 

-0.6202 

-0.8459 

2 

-0.7175 

-1.0118 

3 

-0.7176 

-1.0129 

4 

-0.7176 

-1.0129 

Kant et al [16] 

-0.7164 

- 


Table 6.8: Validation study for moderately thick plate with p refinement for Mesh 

III 

The observations for this study show that: 

• The displacement show convergence to a value of -0.7176 and total potential 
-1.0129 x 10 -4 for III and higher levels of mesh refinement. 

• For a mesh with III refinements, the convergence is observed for p > 3. 





• The difference in the displacement is because in the reference study a fixed 
mesh was used for analysis and no convergence study was done. 

PROBLEM 4 

Four-layer antisymmetric angle-ply square plate subjected to uniformly distributed 
transverse load is analysed for following boundary conditions: 

1. Simply supported (SS) 

2. Clamped edges (CC) 

The laminae properties are: 

fi = 10; E 2 = 6.89 x 10 3 iV/mm 2 ; 

C?i2 = Giz — 0-6 G 2 z = 0.5 E 2 \ 
vn = *d3 = ^23 = 0.25; 

X = 25.4 mm; = 10; 

Thickness of each layer = |; Ply-orientations [45°/ — 45°/45 0 / — 45°]; 
g 0 = 6.89 x lO~ 3 N/mm 2 . 

The results are tabulated in table 6.9, 6.10. 

The observations from our study show that: 

• Both non-dimensionalised displacement and total potential converges for mesh 
refinements above level III. 

• For a fixed level of refinement III, the convergence with p refinement can be 
seen for p > 3. 

• The displacements from our study are higher than those obtained by Reddy 
and -Miravete [29]. We have used HSDT model while they have used CLPT 
and solved the problem by Levy method with state-space approach. 

The HSDT model leads to a more flexible structure as compared to the CLPT 
model. 
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Boundary Condition 

Mesh Type 

Wo 

11^ x 10~ 4 

S S 

Mesh I 

-1.1062 

-2.6566 

Mesh II 

-1.1059 

-2.6682 

Mesh III 

-1.1068 

-2.6698 

Mesh IV 

-1.1069 

-2.6700 

Reddy and Miravete [29] 

- 1.0000 

- 

C C 

Mesh I 

-0.5562 

-1.1342 

Mesh II 

-0.5779 

-1.1756 

Mesh III 

-0.5799 

-1.1819 

Mesh IV 

-0.5800 

-1.1812 

Reddy and Miravete [29] 

-0.5578 

- 


Table 6.9: Validation study for moderately thick plate with mesh refinement for p 

= 3 


Boundary Condition 

p 

Wo 

lip x 10“ 4 

SS 

1 

-0.8862 

-2.0615 

2 



3 

-1.1068 

-2.6698 

4 

-1.1069 

-2.6700 

Reddy and Miravete [29] 

- 1.0000 

- 

C C 

1 

-0.4632 

-0.9143 

2 

-0.5773 

-1.1719 

3 

-0.5799 

-1.1819 

4 

-0.5801 

-1.1$12 

Reddy and Miravete [29] 

-0.5578 

- 


Table 6.10: Validation study for moderately thick plate with p refinement for Mesh 


III 


















6.1.4 THIN PLATES 


Plates with ratio j above 15 are treated to be thin plates. 

PROBLEM 5 

The problem description is same as in problem 1 with the change that £ = 100. 
Our results to this problem are given in table 6.11, 6.12. 


Mesh Type 

Wo 

lip 

Mesh I 

-0.4336 

-0.0605 

Mesh II 

-0.4352 

-0.0612 

Mesh III 

-0.4353 

-0.0614 

Mesh IV 

-0.4353 

-0.0614 

Kant et al [16] 

-0.4344 



Table 6.11: Validation study for thin plate with mesh refinement for p = 3 


V 

Wo 

Ti p 

1 

-0.0388 

-0.0530 

2 

-0.4255 

-0.0600 

3 

-0.4353 

-0.0614 

4 

-0.4353 

-0.0614 

Kant et al [16] 

-0.4344 

- 


Table 6.12: Validation study for thin plate with p refinement for Mesh III 
From above tables we can see that: 

• The non-dimensionalised displacement converges to -0.4353 and total potential 
to -0.0614 for mesh refinement level > III. 


• For fixed mesh refinement level III, and p refinement the results converge with 
p> 3. 





• The non-dimensional displacement from our study is different from the refer- 
ence values. This is because, the reference study was done with a fixed refined 
mesh and convergence study was not done. 

PROBLEM 6 

Simply supported rectangular laminated plate loaded with uniformly distributed 
transverse load is analysed here. The material properties of T300/5208 Graphite/Epoxy 

( pre-preg ) (vj = 0.7) are given below: 

Ei = 132.5 GPa; E 2 = E Z = 10.8 GPa; 

G n = G IZ = 5.7 GPa; G 23 = 3.4 GPa; 
i/ 12 = i/i 3 = 0.24; v 23 = 0.49; 

X = 229 mm; Y = 127 mm; 

Thickness of each ply = 0.127 mm; Ply Orientations [0/90] s ; 

< 7 o = 6.9 x 10~ 4 MPa 

Results are tabulated in table 6.13, 6.14. 

It can be observed that: 

• The displacement converge to value of -11.5004 with mesh refinement for mesh 

refinement level > III. 

• The convergence is with p >3 for p refinement in the mesh with III levels of 
refinement. 

• The reference results are obtained for constant mesh and with FSDT model. 
Also no convergence study was done. We have used HSDT model, which 
assumes flexiblity over FSDT model. 



Mesh Type 

to 0 * 

lip 

Mesh I 


-23.9584 

Mesh II 

-11.5132 

-24.2721 

Mesh III 

-11.5004 

-24.4141 

Mesh IV 



Reddy et al [31] 

-11.5000 

- 


Table 6.13: Validation study for thin plate with mesh refinement for p = 3 


P 

w 0 * 

lip 

1 

-0.0340 

-0.0712 

to 



3 



4 

-11.5004 

-24.4250 

Reddy et al [31] 

-11.5000 

- 


Table 6.14: Validation study for thin plate with p refinement for Mesh III 

















6.2 VALIDATION OF FAILURE CRITERION 


The Non-dimensionalised centre deflection corresponding to the first-ply failure load 

is given by: 


and the non-dimensionalised first-ply failure load ( FLD ) is given by: 

FLD = (f )(|) 4 

For this validation study, mesh topology shown in fig. 6.2 is used. 

First an antisymmetric angle ply laminated plate is studied for its first-ply failure 
load. In second problem, we have analysed symmetric cross ply laminate for its 
first-ply failure load. 

PROBLEM 1 

An antisymmetric angle ply laminated plate clamped along all edges and loaded 
with uniformly distributed transverse load is analysed for first-ply failure load. 

The properties are same as in problem 6 of previous section with the change in the 
ply orientations as [-45/45/ - 45/45]. The material strengths are given as follows: 

X T = 1515 MPa; = 1697 MPa; 

Y t = Z t = 43.8 MPa; Y c = Z c = 43.8 MPa; 

R = 67.6 MPa; S = T = 86.9 MPa 

Results from our analysis are reported in table 6.15 and 6.16. 

In table 6.15-6.18, Rec I and Rec II denotes failure criterion uses the recovered 
stresses along with Eqn. (5.2) and Eqn. (5.3), respectively. Fem denotes finite 
element stresses along with Eqn. (5.3) 

From the above tables we can see that: 

• For fixed value of p with mesh refinement, the failure load for refinement II is 
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very high. This is due to locking effect exhibited by plate. Since plate is very 
thin and clamped along all edges, this effect is more seviour. 

• Similar observation can be made for p refinement for p= 2. 

• With higher levels of refinement in mesh (IV) or high p values the failure loads 
are near to the reference values. 

• The failure loads with recovered stresses are lower compared to failure load 
obtained by using finite element stresses. For the recovered stresses the fail- 
ure criterion yields the same failure load with and with out transverse shear 

stresses. 

• The failure locations are different for our study. The problem considered is 
symmetric in all respect. Hence, locations can also be expected to be symmet- 
ric. We are getting the locations on the opposite face as that obtained from 
reference results. 

• The reference results are obtained for a fixed mesh with FSDT model using 
Lagrangian isoparametric rectangular elements and no convergence study was 

done. 

PROBLEM 2 

The problem description is same as previous problem except the lamination scheme. 
The lamination scheme for this problem is [0/90] s . 

Results from our analysis are reported in table 6.17 and 6.18. 

The comparison of the results with the reference results yield the same conclusions 
as before. 
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Table 6.15: Validation study for first-ply failure load for [-45/45/-45/45] laminate with h 


refinement for p = 2 



Table 6.16: Validation study for first-ply failure load for [-45/45 /-45/45] laminate withp 
refinement for Mesh II 










































































Mesh 

Type 

Stress 

Type 

FLO 

W Q 

xco 

yco 

Ply 

Number 

Location 


Ret I 



m 


4 

Bottom 

Mesh II 




117.06 


4 

Bottom 





119.62 

1.42 

1 

Top 


Rec I 

26955,81 

-21.92 

115.78 


1 

Top 

Mesh III 

Rec II 

26955.81 

-21.92 



1 

Top 


Fern 

29899.92 

-24.32 

117.06 

0.71 

1 

Top 


Rec I 

19499.95 

-18.41 

115.35 

126.05 

1 

Top 

Mesh IV 

Rec II 

19499.95 

-18.41 

115.35 

m 

1 

Top 




-19.85 

112.79 

126.52 

1 

Top 

Bn 


19050.90 

-19.34 



| 

Top 


Table 6.17: Validation study for first-ply failure load for [0/90] s laminate with h refine- 
ment for p = 2 


V 

Stress 

Type 

FLD 

< 

xco 

yco 

Ply 

Number 

Location 


Rec I 

1.215x10® 

-24.38 

117.06 

2.84 

4 

Bottom 

2 

Rec II 

1.215x10® 

-24.38 

117.06 

2.84 

4 

Bottom 


Fern 

1.738x10® 


119.62 


1 

Top 


Rec I 

19423.48 




1 

Top 

3 

Rec II 

19423.48 




1 1 

Top 


Fern 

17588.19 

-17.80 

119.62 

125.57 

1 

Top 


Rec I 

16823.48 

-16.97 

119.62 

1.42 

1 

Top 

4 

Rec II 

16823,48 

-16.97 

119.62 

1.42 

1 

Top 


Fern 

15676.43 

-15.82 

119.62 

125.57 

1 

Top 

Reddy et al [31] 

Fern 

19050.90 

-19.34 

«2 

«75 

1 

Top 


Table 6.18: Validation study for first-ply failure load for [0/90], laminate with p refine- 
ment for Mesh II 


%-7t 'K'CIVTSP? 
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Figure 6.4: Mesh showing the failure points 

In fig. 6.4, the gauss points in element J are shown. The failure indices are calculated 
at these Gauss points, as discussed in the previous chapter. The point 1* and 2* 
denotes the failure locations obtaied by Reddy et al [31] for [-45/45/ — 45/45] and 
[0/90] * laminates, respectively. And 1* and 2* shows the failure locations from our 
study for [-45/45/ - 45/45] and [0/90], laminates, respectively. 
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6.3 OPTIMIZATION 


Following the validation tests of the previous sections, let us now focus on the prob- 
lem of shape optimization. Here, through a few typical examples, we will demon- 
strate the effect of control of discretisation error on the final “optimal” shape of the 
cut-out. 

Note: In each of the problem we start with initial profile with a — b — Yf 8. 
PROBLEM 1: 

A rectangular symmetric laminated plate clamped along one of the smaller edge and 
is loaded transversely with uniformly distributed load with a cut-out is analysed for 
shape optimization. The lamination scheme is [0/90] s . 

Thickness of each layer = 0.127 mm. 
f = 10 Y=f 
q 0 = 2.0N/mm 2 

All laminae are made of T300/5208 Graphite/Epoxy (pre-preg) and their properties 
are given in previous problems. 

The optimum cut-out shape results are tabulated in table 6.19. 



a 

b 

weight 

max. index 

error (in %) 

x coord 

y coord 

Initial 

0.3175 

0.3175 

12.5865 

0.4095 

3.6582 

2.2109 

1.2165 

Intermediate 

0.5799 

0.4109 

12.1544 

0.7709 

4.8672 

1.9517 

1.3472 

Intermediate 

0.6015 

0.4119 

12.1245 

0.7899 

4.3852 

1.9197 

1.2083 

Optimal 

0.6057 

0.4100 

12.1229 

0.7959 

4.8293 

1.9157 

1.2077 


Table 6.19: Optimal shape for cut-out in [0/90] s laminate with out adaptive refine- 

i 

ment and p =2 


From these results we observe that: 

• The errors are below 5% for all the feasible shapes. Hence, adaptive refinement 
is not needed. 
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Mesh I a=b=03175, max, mdex=O.4095, error=3.658% 


Mesh II a=0.5799, b=0.4109, max. index=0.7709, error=4.867 % 



Mesh III »=0.6057, b=0.41, max. index=0.7959, error=4.8293% 



Mesh IV a=Q.6015, b =0.4 119, max. index=0.7899, error=4385% 


Figure 6.5: Cut-out shapes during optimization 


PROBLEM 2: 


The problem description is same as above. The lamination scheme is changed to 
[45/ - 45],. <?o = 0.22 iV/mm 2 


The results are tabulated in table 6.20. 


From these results we observe that: 

• The errors are above 10% in most of the feasible shapes. 

• The final optimum shape is the upper bound on the design variables. 

• The failure occurs at the edges of cut-out boundary and the location changes. 
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E 

b 

weight 

max. index 

error (in %) 

x coord 

y coord 

Initial 

0.3175 

0.3175 

12.5865 

0.4953 

7.1979 

0.0740 

0.0224 

Intermediate 

0.8598 

0.8808 

10.5236 

0.3461 

11.8820 

2.1983 

0.4315 

Intermediate 

1.0122 j 

0.8960 


0.3939 

11.8937 

2.1435 

0.4188 

Optimal 

1.0160 

1.0160 

09.6602 

0.7712 

15.6503 

2.4090 

2.3167 


Table 6.20: Optimal shape for cut-out in [45/ — 45] s laminate with out adaptive 

refinement and p =2 


• The failure is dominated by a yy and r xy . This is in accordance with experi- 
mental study reported by Herakovich [27]. 

For optimum shape obtained, now we will see the effect of adaptive refinement on 
final optimal solution to reduce the error. 

The sequence of adaptive mesh refinement is shown in fig 6.6. The errors for these 
refined meshes are given table 6.21 


Refinement level 

error (in %) 

max. index 

First 



Second 

15.0825 

- 

Third 

12.4125 

- 

Fourth 

11.0823 

- 

Fifth 

10.4591 

1.0089 


Table 6.21: Effect of adaptivity on final optimal shape in [45/ — 45] s laminate 


• Thus, due to adaptive refinement the errors are reduced from 15.6503% to 
10.4591% but the failure index is increased from 0.7712 to 1.0089. 
















Mesh I 


Mesh II 




Mesh III 


Mesh IV 



Mesh V Mesh VI 

Figure 6.6: Adaptive refinement in optimum shape 
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Now we will see the effect of adaptivity when employed through out the optimization 

process. 

The mesh is refined adaptively till either the tolerance in error is reduced below 
7% or for five levels of refinement. The comparison of results with and without 

adaptivity is shown in table 6.22 


Adaptivity 

a 

b 

weight 

max. index 

error (in %) 

No 



09.6602 

0.7712 

15.6503 

Yes 

0.5248 

0.3175 

12.3796 

0.7978 

6.7000 


Table 6.22: Effect of adaptivity on optimal shape in [45/ — 45] s laminate subjected 

to transverse loading 

Thus, it is seen that: 

• The final optimal shape obtained with and without adaptivity are totally dif- 
ferent. 

• The error for the final optimal shape with adaptive procedure is well below 

the tolerance mentioned. 

• The failure occurs in the cut-out boundary but the location changes. 
PROBLEM 3: 

A rectangular symmetrically laminated composite plate clamped along a smaller 
edge and is loaded combinedly with uniformly distributed load, in-plane tensile 
loading and shear load. 

The lamination scheme and laminae properties are as in problem 2. 


q 0 = 0.15 iST/mm 2 ; N x = 20 iV/mm; N xy = 0.35 N/mm 















Adapt ivity 

♦ i | 1 ) 

weight 

max. index 

error (in %) 

No 

1.0155 j 0.8753 

10.1103 

0.7997 

8.3792 

Yes 

1.0157 ! 0.8258 

10.2681 

0.7999 

6.4297 


Table (5.23: Effect of adaptivity an optimal shape in [45/ — 45] s laminate subjected 

to combined loading 

Here, we see t hat : 

• The final optimal shapes obtained with and without adaptivity are near to 

each other. 

• The error in optimum shape obtained with adaptivity is below the specified 

tolerance. 


The failure occurs in the cut-out boundary but the location changes. 



Chapter 7 


CONCLUSIONS AND FUTURE 
SCOPE 

7.1 CONCLUSIONS 

In this study an attempt has been made to estimate and control the discretisation 
error, in the analysis of laminated plates with cut-outs. The effect of errors on the 
evaluated constraint quantities on the final design have also been analysed. From 

this study, it can be further concluded that: 

1. For thick plates, the model used needs moderately refined mesh and approxi 

mation order is sufficient. 

2. For moderately thick plates, proper refinement and approximation orders 

should be chosen for better convergence. 

3. For thin plates, boundary layer effect is severe. Mesh refinements are needed 
in order to eliminate the effect of locking. Higher apporoximation orders (above 
three) is also needed in order to reduce the locking effect. Thus, for a thin plate 
greater care should be taken (with respect to the approximation) in order to 

ensure reliability of the solution quantities. 

4. The recovered stresses seem to be reasonably accurate. A detailed analysis 
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7.2 FUTURE SCOPE 


1. I lie developed finite element program is modular in nature and can, for a given 
mesh, element order p, laminate properties, and loadings, give the solution. 
However, in order to make this code more effective the following modules need 

to he developed: 

(i) Curved elements are to be mapped with higher order approximation so that 
it automatically changes the mapping order according to the approximation 
order used and mapping errors are checked. 

(ii) In the a-posteriori error estimator, recovered stresses are directly used in 
the further analyses. Instead, a more smoothened stress field can be obtained 
by averaging the stresses. 

(iii) The optimization technique uses a simple algorithm of complex search method. 
For sophisticated analysis a suitable Genetic Algorithm can be used. 

(iv) The code is also capable of handling multiple cut-outs. Another module can 
be developed to observe behaviour of plate with multiple cut-outs. 

2. Presently, the code is working well for transverse and in-plane loadings. In 
future modules for buckling analysis, post-buckling behaviour, eigen value prob 
lem, free vibration problems, etc can be added. 

3. Presently optimization uses ‘first-ply failure criterion’ as a constraint. Buck 
ling load can be included as another constraint. 

4. The code is capable of handling multiple load cases. It is to be properly 
activated. 

5. The formulation can be generalized to curvilinear coordinates in order to analyze 
the curved panels and shells also. 

6. Code can be extended to include geometric non-linearity to handle large deformat 
ions. 
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QnHz Q26-H3 QeeHz QnHz QnHz QnHz 
QnHz QnHz QnHz QnH 7 QnH 7 QnH 7 
Q2\Hz Q22HZ Q2&Hz QnH 7 m QnH 7 QuH 7 
QnHz QzeHz QmHz Q\\H 7 QnH 7 QnH 7 


1 



f Quail i (hr t H x QnsHz QibHz 

Q. 15 //i Q.uHi O45H3 QuHz 
QvJh QttHz Q55H5 045^5 
0. !5 // 3 Q44H3 Q45H5 Q- 14^5 


in which, 




1,2, 3, 4, 5,7 


And 


[ft] = 


[Dm] [Do] 0 
[D C ] T [ft] 0 
0 0 [D,] 


XI 



APPENDIX B 




N x 0 0 ... N 2 0 0 

0 Nr 0 ... 0 N 2 0 

0 0 Nr ... 0 0 N 2 


(1) 



Ni, x 

0 

0 

0 

0 

0 

0 


[B»l = 

0 

Ni 

l 'V 

0 

0 

0 

0 

0 

(2) 


Ni 

Ni 

1 ">X 

0 

0 

0 

0 

0 



m = 


0 0 0 Ni 0 0 0 

0 0 0 0 Ni 0 0 

0 0 0 N i<y N, X 0 0 

0 0 0 0 0 Ni, 0 

*sc 

0 0 0 0 0 0 Ni, v 

0 0 0 0 0 N ity N i>x 


(3) 



0 0 N itx 1 0 0 

0 0 TVi 0 10 

0 0 0 0 0 3 

0 0 0 0 0 0 


0 

0 

0 
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(4) 


{d T } = {6l5\ 


2 > 



(5) 
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